library("easypackages")
packages("tidyverse","fpp3", "tsibble", "feasts","fable", "patchwork")
library("USgas")Tarea pronósticos 05
Flujo de trabajo limpio de pronóstico
- Preparación de los datos (limpieza)
- Gráfica de los datos (visualización)
- Definición del modelo (especificación)
- Entrenamiento del modelo (estimación)
- Revisar el desempeño del modelo (evaluación)
- Producir pronósticos
Preparación de los datos (limpieza)
Cargar datos en R, o incluso omitir los datos vacíos (NA), filtrado de la serie (mediante paquetería Tsibble).
Gráfica de los datos (visualización)
global_economy %>%
filter(Country == "Sweden") %>%
autoplot(GDP) +
ggtitle("PIB de Suecia") + ylab("$US billions")Definición del modelo (especificación)
Describir el modelo. Hay distintos tipos, se usará por ejemplo, un modelo lineal de series de tiempo con TLSM:
TSLM(GDP ~ trend())<TSLM model definition>
Entrenamiento del modelo (estimación)
Una vez se define el modelo, lo que sigue es entrenar al modelo. Lo que significa pasarle los datos para que, estdísticamente, encuentre los parámetros que realizan el mejor ajuste posible.
fit <- global_economy %>%
model(Modelo_tendencia = TSLM(GDP ~ trend()))Warning: 7 errors (1 unique) encountered for Modelo_tendencia
[7] 0 (non-NA) cases
Ahora tenemos un modelo lineal ajustado y el objeto resultante es un mable (model table).
fit# A mable: 263 x 2
# Key: Country [263]
Country Modelo_tendencia
<fct> <model>
1 Afghanistan <TSLM>
2 Albania <TSLM>
3 Algeria <TSLM>
4 American Samoa <TSLM>
5 Andorra <TSLM>
6 Angola <TSLM>
7 Antigua and Barbuda <TSLM>
8 Arab World <TSLM>
9 Argentina <TSLM>
10 Armenia <TSLM>
# … with 253 more rows
Revisar el desempeño del modelo (evaluación)
Ver como se ajusta el modelo a los datos y compararlo con otros modelos.
Producir pronósticos
Una vez tenemos un modelo ajustado, podemos hacer pronósticos. Usamos forecast(), h = periodo a pronosticar
fcst <- fit %>% forecast(h = 36) # h = "3 years"
fcst# A fable: 9,468 x 5 [1Y]
# Key: Country, .model [263]
Country .model Year GDP .mean
<fct> <chr> <dbl> <dist> <dbl>
1 Afghanistan Modelo_tendencia 2018 N(1.6e+10, 1.3e+19) 16205101654.
2 Afghanistan Modelo_tendencia 2019 N(1.7e+10, 1.3e+19) 16511878141.
3 Afghanistan Modelo_tendencia 2020 N(1.7e+10, 1.3e+19) 16818654627.
4 Afghanistan Modelo_tendencia 2021 N(1.7e+10, 1.3e+19) 17125431114.
5 Afghanistan Modelo_tendencia 2022 N(1.7e+10, 1.3e+19) 17432207600.
6 Afghanistan Modelo_tendencia 2023 N(1.8e+10, 1.3e+19) 17738984086.
7 Afghanistan Modelo_tendencia 2024 N(1.8e+10, 1.3e+19) 18045760573.
8 Afghanistan Modelo_tendencia 2025 N(1.8e+10, 1.3e+19) 18352537059.
9 Afghanistan Modelo_tendencia 2026 N(1.9e+10, 1.3e+19) 18659313546.
10 Afghanistan Modelo_tendencia 2027 N(1.9e+10, 1.3e+19) 18966090032.
# … with 9,458 more rows
Métodos sencillos de pronóstico
Conocidos como benchmark son básicos pero pueden ser muy utiles.
bricks <- aus_production %>% filter_index("1970" ~ "2004")
bricks# A tsibble: 137 x 7 [1Q]
Quarter Beer Tobacco Bricks Cement Electricity Gas
<qtr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1970 Q1 387 6807 386 1049 12328 12
2 1970 Q2 357 7612 428 1134 14493 18
3 1970 Q3 374 7862 434 1229 15664 23
4 1970 Q4 466 7126 417 1188 13781 20
5 1971 Q1 410 7255 385 1058 13299 19
6 1971 Q2 370 8076 433 1209 15230 23
7 1971 Q3 379 8405 453 1199 16667 28
8 1971 Q4 487 7974 436 1253 14484 24
9 1972 Q1 419 6500 399 1070 13838 24
10 1972 Q2 378 7119 461 1282 15919 34
# … with 127 more rows
bricks %>% autoplot(Bricks)Método del promedio (media)
En este pronósitco, las predicciones de todos los valores futuros son la media de los datos históricos.
bricks %>% model(MEAN(Bricks))# A mable: 1 x 1
`MEAN(Bricks)`
<model>
1 <MEAN>
Método ingenuo (Naive method)
Se toma el último valor como el pronóstico para todos los valores futuros. También se les conoce como pronósitcos de caminata aleatoria, ya que son óptimos en una serie que siga una.
bricks %>% model(NAIVE(Bricks),
RW(Bricks)) # hace exactamente lo mismo que NAIVE()# A mable: 1 x 2
`NAIVE(Bricks)` `RW(Bricks)`
<model> <model>
1 <NAIVE> <NAIVE>
Método ingenuo estacional (seasonal Naive)
Se le agrega un componente a Naive para lidiar con datos altamente estacionales.
bricks %>% model(SNAIVE(Bricks))# A mable: 1 x 1
`SNAIVE(Bricks)`
<model>
1 <SNAIVE>
vic_elec %>%
autoplot(Demand) %>%
plotly::ggplotly()vic_elec %>%
model(SNAIVE(Demand)) %>%
forecast(h = "1 day") %>%
autoplot(vic_elec %>% filter_index("2014-12-30" ~ .))Warning: There was 1 warning in `filter()`.
ℹ In argument: `time_in(Time, ...)`.
Caused by warning:
! Argument 'roll' is deprecated. Deprecated in version '1.8.4'.
Método de la deriva (drift)
Variación de Naive que permite que aumente o disminuya el pronótico en el tiempo. El aumento del cambio es el cambio promedio en los datos históricos.
bricks %>% model(
RW(Bricks ~ drift())
)# A mable: 1 x 1
`RW(Bricks ~ drift())`
<model>
1 <RW w/ drift>
Ejemplos
# Set training data from 1992 to 2006
train <- aus_production %>% filter_index("1992 Q1" ~ "2006 Q4")
# Fit the models
beer_fit <- train %>%
model(
Mean = MEAN(Beer),
`Naïve` = NAIVE(Beer),
`Seasonal naïve` = SNAIVE(Beer),
Drift = RW(Beer ~ drift())
)
# Generate forecasts for 14 quarters
beer_fc <- beer_fit %>% forecast(h= 14)
# Plot forecasts against actual values
beer_fc %>%
autoplot(filter_index(aus_production, "1992 Q1" ~ .), level = NULL) +
ggtitle("Forecasts for quarterly beer production") +
xlab("Year") + ylab("Megalitres") +
guides(colour=guide_legend(title="Forecast")) +
geom_vline(xintercept = as.Date("2007-01-01"), color = "firebrick",
linetype = "dashed") +
annotate("label", x = c(as.Date("2003-01-01"),as.Date("2009-01-01")),
y = 550, label = c("Train set", "Test set"),
color = c("black","blue"))beer_fc %>%
filter(.model == "Seasonal naïve") %>%
autoplot(filter_index(aus_production, "1992 Q1" ~ .)) +
ggtitle("Seasonal naive forecast for quarterly beer production") +
xlab("Quarter") + ylab("Megalitres") +
geom_vline(xintercept = as.Date("2007-01-01"), color = "firebrick", linetype = "dashed") +
annotate("label", x = c(as.Date("2003-01-01"), as.Date("2009-01-01")), y = 550, label = c("Train set", "Test set"), color = c("black", "blue"))Otro ejemplo:
gafa_stock %>% distinct(Symbol)# A tibble: 4 × 1
Symbol
<chr>
1 AAPL
2 AMZN
3 FB
4 GOOG
# Re-index based on trading days
google_stock <- gafa_stock %>%
filter(Symbol == "GOOG") %>%
mutate(day = row_number()) %>%
update_tsibble(index = day, regular = TRUE)
# Filter the year of interest
google_2015 <- google_stock %>% filter(year(Date) == 2015)
# Fit the models
google_fit <- google_2015 %>%
model(
Mean = MEAN(Close),
`Naïve` = NAIVE(Close),
Drift = NAIVE(Close ~ drift()),
SNAIVE = SNAIVE(Close)
)Warning: 1 error encountered for SNAIVE
[1] Non-seasonal model specification provided, use RW() or provide a different lag specification.
google_fit# A mable: 1 x 5
# Key: Symbol [1]
Symbol Mean Naïve Drift SNAIVE
<chr> <model> <model> <model> <model>
1 GOOG <MEAN> <NAIVE> <RW w/ drift> <NULL model>
# Produce forecasts for the 19 trading days in January 2015
google_fc <- google_fit %>% forecast(h = 19)
google_fc# A fable: 76 x 5 [1]
# Key: Symbol, .model [4]
Symbol .model day Close .mean
<chr> <chr> <dbl> <dist> <dbl>
1 GOOG Mean 505 N(602, 6766) 602.
2 GOOG Mean 506 N(602, 6766) 602.
3 GOOG Mean 507 N(602, 6766) 602.
4 GOOG Mean 508 N(602, 6766) 602.
5 GOOG Mean 509 N(602, 6766) 602.
6 GOOG Mean 510 N(602, 6766) 602.
7 GOOG Mean 511 N(602, 6766) 602.
8 GOOG Mean 512 N(602, 6766) 602.
9 GOOG Mean 513 N(602, 6766) 602.
10 GOOG Mean 514 N(602, 6766) 602.
# … with 66 more rows
# A better way using a tsibble to determine the forecast horizons
google_jan_2016 <- google_stock %>%
filter(yearmonth(Date) == yearmonth("2016 Jan"))
google_fc <- google_fit %>% forecast(google_jan_2016)
# Plot the forecasts
google_fc %>%
autoplot(google_2015, level = NULL) +
autolayer(google_jan_2016, Close, color='black') +
ggtitle("Google stock (daily ending 31 Dec 2015)") +
xlab("Day") + ylab("Closing Price (US$)") +
guides(colour=guide_legend(title="Forecast"))Warning: Removed 19 rows containing missing values (`()`).
p1 <- google_fc %>%
filter(.model == "Drift") %>%
autoplot(google_2015) +
autolayer(google_jan_2016, Close, color='black') +
ggtitle("Google stock (daily ending 31 Dec 2015)") +
xlab("Day") + ylab("Closing Price (US$)") +
guides(colour=guide_legend(title="Forecast"))
p2 <- google_fc %>%
filter(.model == "Naïve") %>%
autoplot(google_2015) +
autolayer(google_jan_2016, Close, color='black') +
ggtitle("Google stock (daily ending 31 Dec 2015)") +
xlab("Day") + ylab("Closing Price (US$)") +
guides(colour=guide_legend(title="Forecast"))
p1 / p2Valores ajustados (fitted) y residuales
- Valores reales: Datos históricos de la serie de tiempo.
- Valores ajustados: Valores (fitted) que tratan de pronosticar los valores reales
- Residuales: Valores que el modelo no logró capturar. la diferencia entre valores reales y ajustados.
augment() obtiene los valores ajustados y residuales de un modelo:
augment(beer_fit)# A tsibble: 240 x 6 [1Q]
# Key: .model [4]
.model Quarter Beer .fitted .resid .innov
<chr> <qtr> <dbl> <dbl> <dbl> <dbl>
1 Mean 1992 Q1 443 436. 6.55 6.55
2 Mean 1992 Q2 410 436. -26.4 -26.4
3 Mean 1992 Q3 420 436. -16.4 -16.4
4 Mean 1992 Q4 532 436. 95.6 95.6
5 Mean 1993 Q1 433 436. -3.45 -3.45
6 Mean 1993 Q2 421 436. -15.4 -15.4
7 Mean 1993 Q3 410 436. -26.4 -26.4
8 Mean 1993 Q4 512 436. 75.6 75.6
9 Mean 1994 Q1 449 436. 12.6 12.6
10 Mean 1994 Q2 381 436. -55.4 -55.4
# … with 230 more rows
Diagnóstico de residuales
Es importante analizar los residuales de un modelo. Si los residuales presentan patrones, es signo de que el modelo se puede mejorar.
Los residuales de un buen modelo deben de:
- No estar autocorrelacionadas. Si se detectan correlaciones entre residuos,todavía hay información útil que se debe modelar.
- La media de los residuos es cero. Si es distinta a cero, el pronóstico está sesgado.
Opcionales: - Los residuos tienen una varianza constante. - Los residuos se distribuyen de manera normal.
Transformaciones Box-Cox pueden ayudar a mejorar los residuales
google_2015 %>% autoplot(Close) +
xlab("Day") + ylab("Closing Price (US$)") +
ggtitle("Google Stock in 2015")aug <- google_2015 %>%
model(NAIVE(Close)) %>%
augment()
aug# A tsibble: 252 x 7 [1]
# Key: Symbol, .model [1]
Symbol .model day Close .fitted .resid .innov
<chr> <chr> <int> <dbl> <dbl> <dbl> <dbl>
1 GOOG NAIVE(Close) 253 522. NA NA NA
2 GOOG NAIVE(Close) 254 511. 522. -10.9 -10.9
3 GOOG NAIVE(Close) 255 499. 511. -11.8 -11.8
4 GOOG NAIVE(Close) 256 498. 499. -0.855 -0.855
5 GOOG NAIVE(Close) 257 500. 498. 1.57 1.57
6 GOOG NAIVE(Close) 258 493. 500. -6.47 -6.47
7 GOOG NAIVE(Close) 259 490. 493. -3.60 -3.60
8 GOOG NAIVE(Close) 260 493. 490. 3.61 3.61
9 GOOG NAIVE(Close) 261 498. 493. 4.66 4.66
10 GOOG NAIVE(Close) 262 499. 498. 0.915 0.915
# … with 242 more rows
# aug %>%
# features(.resid, mean(.,na.rm = TRUE))
aug %>% pull(.resid) %>% mean(na.rm = TRUE) [1] 0.9439931
aug %>% autoplot(.resid) + xlab("Día") + ylab("") +
ggtitle("Residuales del método naïve")Warning: Removed 1 row containing missing values (`geom_line()`).
Parece que la media es casi cero y la variación es constante, excepto un outlier.
aug %>%
ggplot(aes(x = .resid)) +
geom_histogram() +
ggtitle("Histograma de los residuales")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
aug %>% ACF(.resid)# A tsibble: 23 x 4 [1]
# Key: Symbol, .model [1]
Symbol .model lag acf
<chr> <chr> <cf_lag> <dbl>
1 GOOG NAIVE(Close) 1 0.0976
2 GOOG NAIVE(Close) 2 -0.0726
3 GOOG NAIVE(Close) 3 -0.0748
4 GOOG NAIVE(Close) 4 -0.0433
5 GOOG NAIVE(Close) 5 -0.0398
6 GOOG NAIVE(Close) 6 0.0206
7 GOOG NAIVE(Close) 7 -0.0652
8 GOOG NAIVE(Close) 8 -0.0291
9 GOOG NAIVE(Close) 9 -0.0379
10 GOOG NAIVE(Close) 10 -0.00687
# … with 13 more rows
aug %>% ACF(.resid) %>% autoplot() + ggtitle("ACF of residuals")Se puede ver que los residuos no están autocorrelacionados (ACF), por lo tanto este método (NAIVE) es bueno para esta serie de tiempo, y los pronósticos derivados de este método pueden ser muy buenos.
Todo en una sola gráfica:
google_2015 %>%
model(NAIVE(Close)) %>%
gg_tsresiduals()Warning: Removed 1 row containing missing values (`geom_line()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
aus_production %>%
filter_index("1992" ~ .) %>%
model(SNAIVE(Beer)) %>%
gg_tsresiduals()Warning: Removed 4 rows containing missing values (`geom_line()`).
Warning: Removed 4 rows containing missing values (`geom_point()`).
Warning: Removed 4 rows containing non-finite values (`stat_bin()`).
aus_production %>%
filter_index("1992" ~ .) %>%
gg_tsdisplay(Beer)Veamos con el método de Media:
google_2015 %>%
model(MEAN(Close)) %>%
augment() %>%
pull(.resid) %>%
mean(na.rm = TRUE)[1] 2.710222e-14
google_2015 %>%
model(MEAN(Close)) %>%
gg_tsresiduals() +
ggtitle("Diagnóstico de residuales para el modelo de Media")Como podemos ver, el método de la media no es óptimo para esta serie de acciones, siendo que la media es 0, al graficar los residuales vemos que muestra un patón (igual al de los datos menos la media), existen claramente autocorrelaciones de una caminata aleatoria y no hay ninguna distribución normal. estos tres factores lo hacen un método nada efectivo para acciones.
Tests de Portmanteau de autocorrelación
Tests para analizar de una forma más formal la presencia o ausencia de autocorrelación en los residuales. Nos ayuda a ver si las primeras \(h\) autocorrelaciones son significativamente distintas de cero o no.
Test de Box-Pierce
Sumatoria de \(r^2_k\) donde \(h\) es el rezago máximo a considerar y \(T\) es la cantidad de observaciones en la muestra. Si \(r^2_k\) es pequeña, entonces el resultado será pequeño. Se recomienda usar \(h = 10\) para datos no estacionales, \(h = 2m\) para datos estaionales (m es el periodo estacional) o como máximo \(h = T/5\).
Hay que tomar en cuenta, que la prueba no es tan buena cuando \(h\) es grande (mayor a \(T/5\)).
Test de Ljung-Box
Generalmente más preciso que Box-Pierce, de igual forma, valores grandes en el resultado indica que no es ruido blanco y los residuales sí están autocorrelacionados.
Hipótesis nula
La hipótesis nula dice que la serie en cuestión no está autocorrelacionada, o sea, es riudo blanco. Para unos residuales, que la hipótesis nula sea cierta es bueno, ya que significa que no están correlacionados los residuales.
¿Cómo saber si es cierta o no?
Siendo \(\alpha\) el nivel de significancia o nivel máximo de error dispuestos a aceptar (usaremos 5%) y p-value un avlor resultante de las pruebas:
Si p-value es menor a \(\alpha\), entonces rechazamos la hiótesis nula, de lo contrario, la aceptamos.
O sea:
- \(p-value > \alpha\), H0 es cierta, residuales no están autocorrelacionados.
- \(p-value < \alpha\), H0 es falsa, residuales están autocorrelacionados.
Ejemplos:
# lag=h and fitdf=K
aug %>% features(.resid, box_pierce, lag=10, dof=0)# A tibble: 1 × 4
Symbol .model bp_stat bp_pvalue
<chr> <chr> <dbl> <dbl>
1 GOOG NAIVE(Close) 7.74 0.654
aug %>% features(.resid, ljung_box, lag=10, dof=0)# A tibble: 1 × 4
Symbol .model lb_stat lb_pvalue
<chr> <chr> <dbl> <dbl>
1 GOOG NAIVE(Close) 7.91 0.637
Como podemos ver para el método Naive, el p-value es mayor al 5%, por lo tanto sus residuos no están autocrrelacionados, son ruido blanco.
google_2015 %>%
model(MEAN(Close)) %>% augment() %>%
features(.resid, box_pierce, lag = 10, dof = 0)# A tibble: 1 × 4
Symbol .model bp_stat bp_pvalue
<chr> <chr> <dbl> <dbl>
1 GOOG MEAN(Close) 2024. 0
google_2015 %>%
model(MEAN(Close)) %>% augment() %>%
features(.resid, ljung_box, lag = 10, dof = 0)# A tibble: 1 × 4
Symbol .model lb_stat lb_pvalue
<chr> <chr> <dbl> <dbl>
1 GOOG MEAN(Close) 2083. 0
Para el caso de Google, no se logra distinguir los residuales del ruido blanco.
Intervalos de predicción
\(c\) es el porcentaje de cobertura de probabilidad. Normalmente se usará 80% y 95%.
Al tener incertidumbre en el ponóstico, entre maoyr sea el intervalo, menos preciso será nuestro pronóstico.
Intervalos de predicción de un paso
Cunado ser realizan pronósitocs de un paso, la desviación estándar del pronóstico es prácticamente la misma que la desviación estándar de los residuos.
Intervalos de predicción de paso múltiple
Entre más aumenta el horizonte de pronóstico, mayor será el intervalo de pronóstico, Esto es, \(\simga_h\) incrementa con \(h\), entonces necesitamos estimaciones de \(\sigma_h\).
Métodos de referencia
Usando fable:
google_2015 %>%
model(NAIVE(Close)) %>%
forecast(h = 10) %>%
hilo()# A tsibble: 10 x 7 [1]
# Key: Symbol, .model [1]
Symbol .model day Close .mean `80%`
<chr> <chr> <dbl> <dist> <dbl> <hilo>
1 GOOG NAIVE(Close) 505 N(759, 125) 759. [744.5400, 773.2200]80
2 GOOG NAIVE(Close) 506 N(759, 250) 759. [738.6001, 779.1599]80
3 GOOG NAIVE(Close) 507 N(759, 376) 759. [734.0423, 783.7177]80
4 GOOG NAIVE(Close) 508 N(759, 501) 759. [730.1999, 787.5601]80
5 GOOG NAIVE(Close) 509 N(759, 626) 759. [726.8147, 790.9453]80
6 GOOG NAIVE(Close) 510 N(759, 751) 759. [723.7543, 794.0058]80
7 GOOG NAIVE(Close) 511 N(759, 876) 759. [720.9399, 796.8202]80
8 GOOG NAIVE(Close) 512 N(759, 1002) 759. [718.3203, 799.4397]80
9 GOOG NAIVE(Close) 513 N(759, 1127) 759. [715.8599, 801.9001]80
10 GOOG NAIVE(Close) 514 N(759, 1252) 759. [713.5329, 804.2272]80
# … with 1 more variable: `95%` <hilo>
google_2015 %>%
model(NAIVE(Close)) %>%
forecast(h = 10) %>%
autoplot(google_2015)Intervalos de predicción con residuales bootstrap
Cuando no es razonable asumir normlidad en los residuos, podemos aplicarles bootstraping, ya que solo asume la no autocorrelación.
Se obteienen muchos escenarios futuros posibles. Para ver algunos, se usa generate:
fit <- google_2015 %>%
model(NAIVE(Close))
sim <- fit %>% generate(h = 30, times = 5, bootstrap = TRUE, seed = 123) # only works with dev version
sim# A tsibble: 150 x 5 [1]
# Key: Symbol, .model, .rep [5]
Symbol .model day .rep .sim
<chr> <chr> <dbl> <chr> <dbl>
1 GOOG NAIVE(Close) 505 1 744.
2 GOOG NAIVE(Close) 506 1 747.
3 GOOG NAIVE(Close) 507 1 733.
4 GOOG NAIVE(Close) 508 1 737.
5 GOOG NAIVE(Close) 509 1 739.
6 GOOG NAIVE(Close) 510 1 733.
7 GOOG NAIVE(Close) 511 1 728.
8 GOOG NAIVE(Close) 512 1 730.
9 GOOG NAIVE(Close) 513 1 721.
10 GOOG NAIVE(Close) 514 1 713.
# … with 140 more rows
Generamos 5 escenarios futuros posibles para los siguientes 30 días de trading.
google_2015 %>%
ggplot(aes(x = day)) +
geom_line(aes(y = Close), size = 1) +
geom_line(aes(y = .sim, colour = as.factor(.rep)), data = sim, size = 1) +
ggtitle("Google closing stock price") +
guides(col = FALSE)Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
of ggplot2 3.3.4.
Con esto, podemos obtener intervalos de predicción, al calcular los percentiles de los escenarios futuros. El resultado es el ntervalo de predicción bootstrapped. Esto se puede lograr con forecast:
fc <- fit %>% forecast(h = 30, bootstrap = TRUE)
fc# A fable: 30 x 5 [1]
# Key: Symbol, .model [1]
Symbol .model day Close .mean
<chr> <chr> <dbl> <dist> <dbl>
1 GOOG NAIVE(Close) 505 sample[5000] 759.
2 GOOG NAIVE(Close) 506 sample[5000] 759.
3 GOOG NAIVE(Close) 507 sample[5000] 759.
4 GOOG NAIVE(Close) 508 sample[5000] 759.
5 GOOG NAIVE(Close) 509 sample[5000] 759.
6 GOOG NAIVE(Close) 510 sample[5000] 759.
7 GOOG NAIVE(Close) 511 sample[5000] 759.
8 GOOG NAIVE(Close) 512 sample[5000] 759.
9 GOOG NAIVE(Close) 513 sample[5000] 759.
10 GOOG NAIVE(Close) 514 sample[5000] 759.
# … with 20 more rows
fc %>% autoplot(google_2015) +
ggtitle("Google closing stock price")fc <- fit %>% forecast(h = 30)
fc# A fable: 30 x 5 [1]
# Key: Symbol, .model [1]
Symbol .model day Close .mean
<chr> <chr> <dbl> <dist> <dbl>
1 GOOG NAIVE(Close) 505 N(759, 125) 759.
2 GOOG NAIVE(Close) 506 N(759, 250) 759.
3 GOOG NAIVE(Close) 507 N(759, 376) 759.
4 GOOG NAIVE(Close) 508 N(759, 501) 759.
5 GOOG NAIVE(Close) 509 N(759, 626) 759.
6 GOOG NAIVE(Close) 510 N(759, 751) 759.
7 GOOG NAIVE(Close) 511 N(759, 876) 759.
8 GOOG NAIVE(Close) 512 N(759, 1002) 759.
9 GOOG NAIVE(Close) 513 N(759, 1127) 759.
10 GOOG NAIVE(Close) 514 N(759, 1252) 759.
# … with 20 more rows
fc %>% autoplot(google_2015) +
ggtitle("Google closing stock price")Pronósticos con transformaciones
Si se le hace un pronóstico a una serie anteriormente transformada (Box-cox por ejemplo), al terminar el pronóstico, se tiene que hacer una transformación inversa para ver la versión original del pronóstico.
Transformación inversa de Box-Cox:
\(y_t=\)
\({exp(w_t), \text{ si }\lambda = 0}.\)
\((\lambda w_t^\lambda+1)^{1/\lambda}, \text{ en otro caso}.\)
La paquetería fable realiza la transofrmación ivnersa en automático, cuando se especifica en la definición del modelo.
Ajustes por sesgo
Un problema al realizar transformaciones matemáticas, como Box–Cox es que las estimaciones puntuales re transformadas ya no presentan la media de la distribución de predicción, sino que la mediana. Esto puede llegar a ser problema dependiendo del contexto, pero normalmente no causa mucha diferencia.
Ajustar por sesgo es la transformación inversa de la media, entre más grande sea la varianza, mayor será la diferencia entre la media y la mediana.
eggs <- as_tsibble(fma::eggs)Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
eggs %>%
model(RW(log(value) ~ drift())) %>%
forecast(h=50) %>%
autoplot(eggs, level = 80,
point_forecast = lst(mean, median))Warning in ggplot2::geom_point(mapping = mapping, data =
dplyr::semi_join(object, : Ignoring unknown aesthetics: linetype
Por lo tanto, para tener un pronóstico sesgado, se usa point_forecast = lst(median)
Pronósticos con descomposición
La descomposición de series de tiempo puede ser útil para producir pronósticos.
El pronóstico se realiza en dos pasos: - Un pronóstico para el componente estacional. - Un pronósitco separado para la serie desestacionalizada.
En realidad, SNAIVE es el pronóstido de solo el componente estacional.
(us_retail_employment <- fpp3::us_employment %>%
filter(year(Month) >= 1990, Title == "Retail Trade"))# A tsibble: 357 x 4 [1M]
# Key: Series_ID [1]
Month Series_ID Title Employed
<mth> <chr> <chr> <dbl>
1 1990 Jan CEU4200000001 Retail Trade 13256.
2 1990 Feb CEU4200000001 Retail Trade 12966.
3 1990 Mar CEU4200000001 Retail Trade 12938.
4 1990 Apr CEU4200000001 Retail Trade 13012.
5 1990 May CEU4200000001 Retail Trade 13108.
6 1990 Jun CEU4200000001 Retail Trade 13183.
7 1990 Jul CEU4200000001 Retail Trade 13170.
8 1990 Aug CEU4200000001 Retail Trade 13160.
9 1990 Sep CEU4200000001 Retail Trade 13113.
10 1990 Oct CEU4200000001 Retail Trade 13185.
# … with 347 more rows
us_retail_employment %>%
autoplot(Employed)dcmp <- us_retail_employment %>%
model(STL(Employed ~ trend(window = 7), robust=TRUE)) %>%
components() %>%
select(-.model)
dcmp# A tsibble: 357 x 7 [1M]
# Key: Series_ID [1]
Series_ID Month Employed trend season_year remainder season_adjust
<chr> <mth> <dbl> <dbl> <dbl> <dbl> <dbl>
1 CEU4200000001 1990 Jan 13256. 13326. -72.6 2.85 13328.
2 CEU4200000001 1990 Feb 12966. 13297. -273. -58.1 13239.
3 CEU4200000001 1990 Mar 12938. 13270. -295. -36.3 13233.
4 CEU4200000001 1990 Apr 13012. 13230. -216. -1.59 13228.
5 CEU4200000001 1990 May 13108. 13223. -119. 4.72 13228.
6 CEU4200000001 1990 Jun 13183. 13212. -30.6 0.978 13213.
7 CEU4200000001 1990 Jul 13170. 13198. -32.8 5.23 13203.
8 CEU4200000001 1990 Aug 13160. 13178. -20.1 1.84 13180.
9 CEU4200000001 1990 Sep 13113. 13149. -45.5 9.76 13159.
10 CEU4200000001 1990 Oct 13185. 13113. 64.2 8.28 13121.
# … with 347 more rows
dcmp %>%
model(NAIVE(season_adjust)) %>%
forecast() %>%
autoplot(dcmp) + ylab("New orders index") +
ggtitle("Pronóstico naïve de los datos desestacionalizados")Le podemos agregar nuevamente la estacionalidad con la función decomposition_model():
us_retail_employment %>%
model(stlf = decomposition_model(
STL(Employed ~ trend(window = 7), robust = TRUE),
NAIVE(season_adjust)
)) %>%
forecast() %>%
autoplot(us_retail_employment)Evaluación del desempeño de los pronósticos
Conjuntos de entrenamiento y prueba
El tamaño de la prueba es, generalmente, del 20% del total de datos disponibles.
- Un modelo que se ajusta muy bien a los datos de entrenamiento no necesariamente produce los mejore spronósticos.
- Al aumentar los parámetros, podemos llegar a un ajuste perfecto del modelo a los datos.
- Puede darse un efecto de sobre ajuste (over-fitting) y esto es tan malo como tener un muy mal ajuste.
Funciones para sementar las series de tiempo
top_n nos permite obtener las n observaciones más extremas.
gafa_stock %>%
group_by(Symbol) %>%
top_n(1, Close)# A tsibble: 4 x 8 [!]
# Key: Symbol [4]
# Groups: Symbol [4]
Symbol Date Open High Low Close Adj_Close Volume
<chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 AAPL 2018-10-03 230. 233. 230. 232. 230. 28654800
2 AMZN 2018-09-04 2026. 2050. 2013 2040. 2040. 5721100
3 FB 2018-07-25 216. 219. 214. 218. 218. 58954200
4 GOOG 2018-07-26 1251 1270. 1249. 1268. 1268. 2405600
Errores de pronóstico
Diferencia entre el valor real ocurrido y el dato pronosticado.
Los errores de pronóstico son distintos de los residuales en:
- Los residuales se calculan con los datos de entrenamiento, mientras que los errores de pronóstico se calculan con los de prueba.
- Los residuos se calculan mediante pronósticos de un paso, los errores pueden ser multi-step.
Errores dependientes de la escala de los datos
Dependen de la escala de la serie de tiempo, por lo que no se puede comparar con series de tiempo con otras unidades.
Los dos más utilizados son:
- MAE, Mean absolute error, el promedio del error.
- RMSE, Root mean squared error, la raíz cuadrada del promedio del error.
MAE es muy utilizado por su facilidad de cómputo e interpretación.
Errores porcentuales
Son un porcentaje y son muy utilizados para comparar el desempeño de pronósticos de distintos conjuntos de datos.
-MAPE, Mean absolute percentage error, el promedio del error.
El problema es que si un dato es 0, este error será infinito, por lo tanto se inventó:
-sMAPE, symmetric MAPE, resuelve el problema del MAPE.
Sin embargo, no es recomendable en la práctica.
Errores escalados
ALternativa para porcentuales para comparar con otras series.
- MASE, Mean absolute scalated error, promedio de los errores escalados.
recent_production <- aus_production %>% filter(year(Quarter) >= 1992)
beer_train <- recent_production %>% filter(year(Quarter) <= 2007)
beer_fit <- beer_train %>%
model(
Mean = MEAN(Beer),
`Naïve` = NAIVE(Beer),
`Seasonal naïve` = SNAIVE(Beer),
Drift = RW(Beer ~ drift())
)
beer_fc <- beer_fit %>%
forecast(h = 10)
beer_fc %>%
autoplot(recent_production, level = NULL) +
xlab("Year") + ylab("Megalitres") +
ggtitle("Forecasts for quarterly beer production") +
guides(colour=guide_legend(title="Forecast"))Usamos accuracy para ver los errores en el ajuste de modelos (mabble) a los datos de entrenamiento:
accuracy(beer_fit)# A tibble: 4 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Mean Training 0 43.6 35.2 -0.937 7.89 2.46 2.60 -0.109
2 Naïve Training 4.76e- 1 65.3 54.7 -0.916 12.2 3.83 3.89 -0.241
3 Seasonal naïve Training -2.13e+ 0 16.8 14.3 -0.554 3.31 1 1 -0.288
4 Drift Training -5.41e-15 65.3 54.8 -1.03 12.2 3.83 3.89 -0.241
Ahora la usamos para ver los errores de pronóstico (fabble) comparado a los datos de prueba:
beer_fc |> accuracy(recent_production)# A tibble: 4 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Drift Test -54.0 64.9 58.9 -13.6 14.6 4.12 3.87 -0.0741
2 Mean Test -13.8 38.4 34.8 -3.97 8.28 2.44 2.29 -0.0691
3 Naïve Test -51.4 62.7 57.4 -13.0 14.2 4.01 3.74 -0.0691
4 Seasonal naïve Test 5.2 14.3 13.4 1.15 3.17 0.937 0.853 0.132
Los errores de ajuste (datos de entrenamiento) sirven para ver que modelos utilizar para pronosticar.
Los errores de pronóstico (datos de prueba) nos da mayor claridad de cuáles modelos parecen producir los mejores pronósticos.
Después de esto se recalculan los modelos utilizando toda la información (sin separar en conjuntos de entrenamiento y prueba) y se producen los pronósticos reales a presentar.
Tarea
Series:
1. Serie de ventas de ropa en el Estado de Victoria, Australia
(clothes_victoria <- filter(aus_retail, State == "Victoria", Industry == "Clothing retailing") |> select(Month, Turnover))# A tsibble: 441 x 2 [1M]
Month Turnover
<mth> <dbl>
1 1982 Apr 93.6
2 1982 May 95.3
3 1982 Jun 85.2
4 1982 Jul 91.6
5 1982 Aug 85.2
6 1982 Sep 89.5
7 1982 Oct 93
8 1982 Nov 108.
9 1982 Dec 148.
10 1983 Jan 81.6
# … with 431 more rows
train <- filter_index(clothes_victoria, "2000 Jan" ~ "2015 Dec")
autoplot(clothes_victoria, Turnover)fit <- model(train,
Mean = MEAN(Turnover),
Naive = NAIVE(Turnover),
Snaive = SNAIVE(Turnover),
Drift = RW(Turnover ~ drift()))
fc <- forecast(fit, h = "3 years")
fc |> autoplot(filter_index(clothes_victoria, "2000 Jan" ~ .), level = NULL) +
ggtitle("Pronóstico de ventas mensuales de ropa en Victoria") +
xlab("Mes") + ylab("MegaWhats/hora") +
guides(color=guide_legend(title="Tipo")) +
geom_vline(xintercept = as.Date("2016-01-01"), color = "firebrick", linetype = "dashed")select(fit, Mean) |> gg_tsresiduals() + ggtitle("Media") La media no es buena idea, no muestra distribución normal del todo, están autocorrelacionados y la gráfica muestra un patrón.
select(fit, Naive) |> gg_tsresiduals() + ggtitle("Naive") Warning: Removed 1 row containing missing values (`geom_line()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
Naive parece ser la mejor opción, tiene distribución normal, su varianza es algo constante y no parece estar muy correlacionado.
select(fit, Snaive) |> gg_tsresiduals() + ggtitle("Snaive") Warning: Removed 12 rows containing missing values (`geom_line()`).
Warning: Removed 12 rows containing missing values (`geom_point()`).
Warning: Removed 12 rows containing non-finite values (`stat_bin()`).
Snaive muestra distribución normal, sin embargo, parece presentar cierta correlación y patrón. A pesar de esto, Snaive sigue siendo un buen modelo al ver el gráfico.
select(fit, Drift) |> gg_tsresiduals() + ggtitle("Drift") Warning: Removed 1 row containing missing values (`geom_line()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
Drift tiene residuales muy similares a Naive, pueden ser buenas bases para aplicar otros modelos, sin embargo, tomando solo en cuenta los benchmark, probablemente Snaive será el que se va a utilizar.
fc <- filter(fc, .model == "Snaive")
fc_bs <- forecast(fit, h ="3 years", bootstrap = TRUE)
fc_bs <- filter(fc_bs, .model == "Snaive")
fc |> autoplot(filter_index(clothes_victoria, "2000 Jan" ~ .)) +
ggtitle("Pronóstico de ventas mensuales de ropa en Victoria") +
xlab("Mes") + ylab("MegaWhats/hora") +
guides(color=guide_legend(title="Tipo")) +
geom_vline(xintercept = as.Date("2016-01-01"), color = "firebrick", linetype = "dashed")fc_bs |> autoplot(filter_index(clothes_victoria, "2000 Jan" ~ .)) +
ggtitle("Pronóstico de ventas mensuales de ropa en Victoria") +
xlab("Mes") + ylab("MegaWhats/hora") +
guides(color=guide_legend(title="Tipo")) +
geom_vline(xintercept = as.Date("2016-01-01"), color = "firebrick", linetype = "dashed")Al no ser datos muy aleatorios o los datos de una acción, y que los horizontes de pronóstico no son muy diferentes, no es necesario usar Bootstrap.
2. Serie de PIB per capita de Australia
(aus_gdp <- filter(global_economy, Code == "AUS") %>%
mutate(GDP_pc = GDP / Population) |> select(Year, GDP_pc))# A tsibble: 58 x 2 [1Y]
Year GDP_pc
<dbl> <dbl>
1 1960 1807.
2 1961 1874.
3 1962 1851.
4 1963 1964.
5 1964 2128.
6 1965 2277.
7 1966 2340.
8 1967 2576.
9 1968 2719.
10 1969 2986.
# … with 48 more rows
train <- filter_index(aus_gdp, "1970" ~ "2010")
autoplot(aus_gdp, GDP_pc)fit <- model(train,
Mean = MEAN(GDP_pc),
Naive = NAIVE(GDP_pc),
Snaive = SNAIVE(GDP_pc),
Drift = RW(GDP_pc ~ drift()))Warning: 1 error encountered for Snaive
[1] Non-seasonal model specification provided, use RW() or provide a different lag specification.
fc <- forecast(fit, h = "7 years")
fc |> autoplot(filter_index(aus_gdp, "1970" ~ .), level = NULL) +
ggtitle("Pronóstico de PIB per Capita de Australia") +
xlab("Año") + ylab("USD") +
guides(color=guide_legend(title="Tipo")) +
geom_vline(xintercept = 2011, color = "firebrick", linetype = "dashed")Warning: Removed 7 rows containing missing values (`()`).
Como podemos ver, Snaive no funciona al no ser estacional (se usan datos anuales)
select(fit, Mean) |> gg_tsresiduals() + ggtitle("Media") La media no es buena idea, no muestra distribución normal, están autocorrelacionados y la gráfica muestra un patrón muy similar a los datos.
select(fit, Naive) |> gg_tsresiduals() + ggtitle("Naive") Warning: Removed 1 row containing missing values (`geom_line()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
Naive parece ser la mejor opción, tiene distribución normal, y a pesar de no tener una varianza muy constante, es la mejor, y no parece estar muy correlacionado.
select(fit, Drift) |> gg_tsresiduals() + ggtitle("Drift") Warning: Removed 1 row containing missing values (`geom_line()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
Drift tiene residuales muy similares a Naive, y tomando en cuenta las gráficas, se decide usar Drift para capturar la tendencia alcista que tiene.
Al no ser datos muy aleatorios o los datos de una acción, y que los horizontes de pronóstico no son muy diferentes, no es necesario usar Bootstrap.